clear all; clear global; clc; close all;

load('../Set_postestimation_use.mat');

tau_A = 0.053;
tau_B = 0.038;

tau_bA = 0.053;
tau_bB = 0.041;

Par.tau_A = tau_A;
Par.tau_B = tau_B;

Par.L_A   = 1;
Par.L_B   = 1;

Par.time_ss  = 2*Par.time_ss;
Par.time_81  = Par.time1;
Par.time_95  = 20/Par.dt;
Par.time_add = (200+Par.time_81*Par.dt)/Par.dt-Par.time_ss+1;

Par.num_r = Par.num_r/2;

horizons = 5:5:50;
tau_bA_grid = 0.35; %linspace(0,0.60,301);

[Par,Step,Trans] = Trade_ctrf(Par);

rev_b   = 1; % Revenue back to consumer
retal_b = 0; % 1 if foreign retaliation

Par.rev_US = rev_b;
Par.rev_FN = rev_b;

% ============================================
%                 SOLUTIONS
% ============================================

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

steps      = Step{1};
Trans_innov = Step{2};
Trans_fail = Step{3};
Trans_exo  = Step{4};

%% Calibrated Economy Path
% =================================================================

Tol_rr = 1*10^-6;
save Tol_rr Tol_rr

%---------------- Path

mu_init = y_init;
Q_init  = y_init;

Par.if_long = 0;

[rr_sqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, Ql_sqn, mu_share, ind_neg, ...
                            prft_fnc_sqn, wage_dom_sqn, mcuts] = ...
                                Path_ctrf(mu_init,Q_init,Par,Step,Trans);
                            
%---------------- Variables

% Profits
% --------------------------------------------------------

prft_fnc_A = prft_fnc_sqn(1:size(prft_fnc_sqn,1)/2,:);
prft_fnc_B = prft_fnc_sqn(size(prft_fnc_sqn,1)/2+1:end,:);

wage_dom_A = wage_dom_sqn(1:size(wage_dom_sqn,1)/2,:);
wage_dom_B = wage_dom_sqn(size(wage_dom_sqn,1)/2+1:end,:);

% Qualities
% --------------------------------------------------------

QA = Q_sqn(1:size(Q_sqn,1)/2,:);
QB = Q_sqn(size(Q_sqn,1)/2+1:end,:);

QA_long = Ql_sqn(1:size(Q_sqn,1)/2,:);
QB_long = Ql_sqn(size(Q_sqn,1)/2+1:end,:);

QwA = sum(wage_dom_A.*QA+flip(1-wage_dom_B).*QA*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet));
QwB = sum(wage_dom_B.*QB+flip(1-wage_dom_A).*QB*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet));

qq_sqn = [sum(QA)./QwA; sum(QB)./QwB];
                
wq_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(1,:).^(-Par.bet);
            (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(2,:).^(-Par.bet)];  

xk_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wq_sqn(1,:).^(1/Par.bet)*Par.L_B;
         (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wq_sqn(2,:).^(1/Par.bet)*Par.L_A];

% Shares
% --------------------------------------------------------

MU      = mu_share(:,1:time_ss+1);
MU_long = mu_share;

McutX  = mcuts(1,1:time_ss);
McutM  = mcuts(2,1:time_ss);

% R&D
% --------------------------------------------------------

X_iA = X_sqn(1:size(X_sqn,1)/2,:);
X_iB = X_sqn(size(X_sqn,1)/2+1:end,:);

V_iA = V_sqn(1:size(V_sqn,1)/2,:);
V_iB = V_sqn(size(V_sqn,1)/2+1:end,:);

X_eA = Xe_sqn(1:size(Xe_sqn,1)/2,:);
X_eB = Xe_sqn(size(Xe_sqn,1)/2+1:end,:);

V_eA = alfa_eA*X_eA.^gama_eA*(1-1/gama_eA);
V_eB = alfa_eB*X_eB.^gama_eB*(1-1/gama_eB);

% Income
% --------------------------------------------------------

II.IA_1_pr  = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet);
II.IA_2_rdi = sum(alfa_A/gama_A*X_iA.^gama_A.*QA);
II.IA_3_rde = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA);
II.IA_4_wd  = pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA);
II.IA_5_wf  = pi_frn_toUS/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB);
II.IA_6_wig = wq_sqn(1,:).*sum(QA);
II.IA_7_trf = rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);

II.IB_1_pr  = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet);
II.IB_2_rdi = sum(alfa_B/gama_B*X_iB.^gama_B.*QB);
II.IB_3_rde = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB);
II.IB_4_wd  = pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB);
II.IB_5_wf  = pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA);
II.IB_6_wig = wq_sqn(2,:).*sum(QB);
II.IB_7_trf = rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA);

IA = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet)-...
    sum(alfa_A/gama_A*X_iA.^gama_A.*QA)-...
    sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA)+...
    pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA)+...
    pi_frn_toUS/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB)+...
    wq_sqn(1,:).*sum(QA)+...
    rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);  
IB = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet)-...
    sum(alfa_B/gama_B*X_iB.^gama_B.*QB)-...
    sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB)+...
    pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB)+...
    pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA)+...
    wq_sqn(2,:).*sum(QB)+...
    rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA); 

g_IA = mean((IA(2:time_81+1)-IA(1:time_81))./IA(1:time_81));
g_IB = mean((IB(2:time_81+1)-IB(1:time_81))./IB(1:time_81)); 

gI_A = (IA(2:end)./IA(1:end-1))-1;
gI_B = (IB(2:end)./IB(1:end-1))-1;


% Welfare
% --------------------------------------------------------

if psy == 1
    
    for ctr_t = 1:length(horizons)
    
        Wlf_o(ctr_t,1) = sum(exp(-ro*(dt:dt:horizons(ctr_t))).*log(IA(time_81+1:time_81+horizons(ctr_t)/dt)));

        Wlf_oB(ctr_t,1) = sum(exp(-ro*(dt:dt:horizons(ctr_t))).*log(IB(time_81+1:time_81+horizons(ctr_t)/dt)));
        
    end  
    
else
    
    for ctr_t = 1:length(horizons)
        
        Wlf_o(ctr_t,1) = (1-psy)^-1*[sum(exp(-ro*(dt:dt:horizons(ctr_t))).*(IA(time_81+1:time_81+horizons(ctr_t)/dt).^(1-psy)))];

        Wlf_oB(ctr_t,1) = (1-psy)^-1*[sum(exp(-ro*(dt:dt:horizons(ctr_t))).*(IB(time_81+1:time_81+horizons(ctr_t)/dt).^(1-psy)))];
        
    end  
   
end
  

%% Policy Change
% =================================================================
%

Tol_rr = 1*10^-6;
save Tol_rr Tol_rr

% ------------- Parfor intialization ------------- 

mub_init = mu_share(:,time_81+1);
Qb_init  = QA(:,time_81+1);

mu_pre = mu_share(:,1:time_81);
rr_pre = rr_sqn(:,1:time_81);
IA_pre = IA(:,1:time_81);
IB_pre = IB(:,1:time_81);
QA_pre = QA(:,1:time_81);
QB_pre = QB(:,1:time_81);
QwA_pre = QwA(:,1:time_81);
QwB_pre = QwB(:,1:time_81);
McutM_pre = McutM(1,1:time_81);
McutX_pre = McutX(1,1:time_81);

QA_post = QA(:,time_81+1:end);
QB_post = QB(:,time_81+1:end);

II_orig = II;

MCUTX = [McutX; zeros(length(tau_bA_grid),size(McutX,2))];
MCUTM = [McutM; zeros(length(tau_bA_grid),size(McutM,2))];

% parpool(24)

for ctr_s = 1:length(tau_bA_grid)
    
    % =============== Policy Change ===================

    Par_b = Par;
    Par_b.num_r = Par.num_r/2;

    Par_b.tau_A = tau_bA_grid(ctr_s);
    Par_b.tau_B = tau_bB;
    
    Par_b.rev_US = rev_b;
    Par_b.rev_FN = rev_b;

    [Par_b,Step_b,Trans_b] = Trade_ctrf(Par_b);

    % ============================================
    %                 SOLUTIONS
    % ============================================

    pi_domb = Par_b.pi_dom;
    pi_frnb = Par_b.pi_frn;
    pi_frnb_toUS = Par_b.pi_frn_toUS;

    %---------------- Path

    [rrb_sqn, Xb_sqn, Vb_sqn, Xeb_sqn, Qb_sqn, Qlb_sqn, mub_share, ind_neg_b, ...
                                prft_fncb_sqn, wage_domb_sqn, mcutsb] = ...
                                Path_ctrf(mub_init,Qb_init,Par_b,Step_b,Trans_b);

    rr_b = [rr_pre rrb_sqn(:,1:end-time_81)];

    %---------------- Variables

    % Profits
    % --------------------------------------------------------

    prft_fncb_A = prft_fncb_sqn(1:size(prft_fncb_sqn,1)/2,:);
    prft_fncb_B = prft_fncb_sqn(size(prft_fncb_sqn,1)/2+1:end,:);

    wage_domb_A = wage_domb_sqn(1:size(wage_domb_sqn,1)/2,:);
    wage_domb_B = wage_domb_sqn(size(wage_domb_sqn,1)/2+1:end,:);

    % Qualities
    % --------------------------------------------------------

    QA_b = Qb_sqn(1:size(Q_sqn,1)/2,:);
    QB_b = Qb_sqn(size(Q_sqn,1)/2+1:end,:);

    QAb_long = Qlb_sqn(1:size(Q_sqn,1)/2,:);
    QBb_long = Qlb_sqn(size(Q_sqn,1)/2+1:end,:);

    QwA_b = sum(wage_domb_A.*QA_b+flip(1-wage_domb_B).*QA_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
    QwB_b = sum(wage_domb_B.*QB_b+flip(1-wage_domb_A).*QB_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

    qqb_sqn = [sum(QA_b)./QwA_b; sum(QB_b)./QwB_b];

    wqb_sqn = [(1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(1,:).^(-Par_b.bet);
                (1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(2,:).^(-Par_b.bet)];  

    xkb_sqn = [(Par_b.pi_frn/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(1,:).^(1/Par_b.bet)*Par.L_B;
                (Par_b.pi_frn_toUS/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(2,:).^(1/Par_b.bet)*Par.L_A];

    QA_orgb = [QA_pre QA_b(:,1:end-time_81)];
    QB_orgb = [QB_pre QB_b(:,1:end-time_81)];

    QwA_orgb = [QwA_pre QwA_b(:,1:end-time_81)];
    QwB_orgb = [QwB_pre QwB_b(:,1:end-time_81)];

    qq_orgb  = [sum(QA_orgb)./QwA_orgb; sum(QB_orgb)./QwB_orgb];

    % Shares
    % --------------------------------------------------------

    MU_b = [mu_pre mub_share(:,1:end-time_81)];

    MCUTX(ctr_s+1,:)  = [McutX_pre mcutsb(1,1:time_ss-time_81)];
    MCUTM(ctr_s+1,:)  = [McutM_pre mcutsb(2,1:time_ss-time_81)];

    % R&D
    % --------------------------------------------------------

    Xb_iA = Xb_sqn(1:size(X_sqn,1)/2,:);
    Xb_iB = Xb_sqn(size(X_sqn,1)/2+1:end,:);

    Xb_eA = Xeb_sqn(1:size(Xe_sqn,1)/2,:);
    Xb_eB = Xeb_sqn(size(Xe_sqn,1)/2+1:end,:);

    % Income 
    % --------------------------------------------------------

    II = II_orig;

    II.IA_bb1_pr  = sum(prft_fncb_A.*QA_b).*wqb_sqn(1,:).^(1-1/bet);
    II.IA_bb2_rdi = sum(alfa_A/gama_A*Xb_iA.^gama_A.*QA_b);
    II.IA_bb3_rde = sum(alfa_eA/gama_eA*Xb_eA.^gama_eA.*QA_b);
    II.IA_bb4_wd  = pi_domb/(1-bet)*L_A*wqb_sqn(1,:).^(1-1/bet).*sum(wage_domb_A.*QA_b);
    II.IA_bb5_wf  = pi_frnb_toUS/(1-bet)*L_A*wqb_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_domb_A).*QB_b);
    II.IA_bb6_wig = wqb_sqn(1,:).*sum(QA_b);
    II.IA_bb7_trf = Par_b.rev_US*Par_b.trf_A*((1+Par_b.trf_A)*(1+kapa)*eta/(1-bet))*wqb_sqn(2,:).*xkb_sqn(2,:).*sum(flip(1-wage_domb_A).*QB_b);

    II.IA_bb2_xi = sum(alfa_A/gama_A*X_iA.^gama_A.*[QA(:,1:time_81) QA_b(:,1:end-time_81)]);
    II.IA_bb3_xe = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*[QA(:,1:time_81) QA_b(:,1:end-time_81)]);

    II.IB_bb1_pr  = sum(prft_fncb_B.*QB_b).*wqb_sqn(2,:).^(1-1/bet);
    II.IB_bb2_rdi = sum(alfa_B/gama_B*Xb_iB.^gama_B.*QB_b);
    II.IB_bb3_rde = sum(alfa_eB/gama_eB*Xb_eB.^gama_eB.*QB_b);
    II.IB_bb4_wd  = pi_domb/(1-bet)*L_B*wqb_sqn(2,:).^(1-1/bet).*sum(wage_domb_B.*QB_b);
    II.IB_bb5_wf  = pi_frnb/(1-bet)*L_B*wqb_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_domb_B).*QA_b);
    II.IB_bb6_wig = wqb_sqn(2,:).*sum(QB_b);
    II.IB_bb7_trf = Par_b.rev_FN*Par_b.trf_B*((1+Par_b.trf_B)*(1+kapa)*eta/(1-bet))*wqb_sqn(1,:).*xkb_sqn(1,:).*sum(flip(1-wage_domb_B).*QA_b);

    II.IB_bb2_xi = sum(alfa_B/gama_B*X_iB.^gama_B.*[QB(:,1:time_81) QB_b(:,1:end-time_81)]);
    II.IB_bb3_xe = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*[QB(:,1:time_81) QB_b(:,1:end-time_81)]);

    IA_bb = II.IA_bb1_pr-II.IA_bb2_rdi-II.IA_bb3_rde+...
            II.IA_bb4_wd+II.IA_bb5_wf+II.IA_bb6_wig+II.IA_bb7_trf;

    IB_bb = II.IB_bb1_pr-II.IB_bb2_rdi-II.IB_bb3_rde+...
            II.IB_bb4_wd+II.IB_bb5_wf+II.IB_bb6_wig+II.IB_bb7_trf;

    IA_b = [IA_pre IA_bb(:,1:end-time_81)];
    IB_b = [IB_pre IB_bb(:,1:end-time_81)];

    % Alternative for Welfare

    QA_bo = [QA_post 0*QA_b(:,end-time_81+1:end)];
    QB_bo = [QB_post 0*QB_b(:,end-time_81+1:end)];

    QwA_bo = sum(wage_domb_A.*QA_bo+flip(1-wage_domb_B).*QA_bo*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
    QwB_bo = sum(wage_domb_B.*QB_bo+flip(1-wage_domb_A).*QB_bo*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

    qqbo_sqn = [sum(QA_bo)./QwA_bo; sum(QB_bo)./QwB_bo];

    wqbo_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qqbo_sqn(1,:).^(-Par.bet);
                (1-Par.bet)*Par.eta^(Par.bet-1)*qqbo_sqn(2,:).^(-Par.bet)];  

    xkbo_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wqbo_sqn(1,:).^(1/Par.bet)*Par.L_B;
                (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wqbo_sqn(2,:).^(1/Par.bet)*Par.L_A];

    II.IA_bo1_pr  = sum(prft_fncb_A.*QA_bo).*wqbo_sqn(1,:).^(1-1/bet);
    II.IA_bo2_rdi = sum(alfa_A/gama_A*Xb_iA.^gama_A.*QA_bo);
    II.IA_bo3_rde = sum(alfa_eA/gama_eA*Xb_eA.^gama_eA.*QA_bo);
    II.IA_bo4_wd  = pi_domb/(1-bet)*L_A*wqbo_sqn(1,:).^(1-1/bet).*sum(wage_domb_A.*QA_bo);
    II.IA_bo5_wf  = pi_frnb_toUS/(1-bet)*L_A*wqbo_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_domb_A).*QB_bo);
    II.IA_bo6_wig = wqbo_sqn(1,:).*sum(QA_bo);
    II.IA_bo7_trf = Par_b.rev_US*Par_b.trf_A*((1+Par_b.trf_A)*(1+kapa)*eta/(1-bet))*wqbo_sqn(2,:).*xkbo_sqn(2,:).*sum(flip(1-wage_domb_A).*QB_bo);

    II.IB_bo1_pr  = sum(prft_fncb_B.*QB_bo).*wqbo_sqn(2,:).^(1-1/bet);
    II.IB_bo2_rdi = sum(alfa_B/gama_B*Xb_iB.^gama_B.*QB_bo);
    II.IB_bo3_rde = sum(alfa_eB/gama_eB*Xb_eB.^gama_eB.*QB_bo);
    II.IB_bo4_wd  = pi_domb/(1-bet)*L_B*wqbo_sqn(2,:).^(1-1/bet).*sum(wage_domb_B.*QB_bo);
    II.IB_bo5_wf  = pi_frnb/(1-bet)*L_B*wqbo_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_domb_B).*QA_bo);
    II.IB_bo6_wig = wqbo_sqn(2,:).*sum(QB_bo);
    II.IB_bo7_trf = Par_b.rev_FN*Par_b.trf_B*((1+Par_b.trf_B)*(1+kapa)*eta/(1-bet))*wqbo_sqn(1,:).*xkbo_sqn(1,:).*sum(flip(1-wage_domb_B).*QA_bo);

    % Welfare
    % --------------------------------------------------------

    [WlfA,WlfB,IAb_set,IBb_set] = postestimation_welfare_trade(II,Par,horizons); 

    WLFA(:,ctr_s) = WlfA(:,1);
    WLFB(:,ctr_s) = WlfB(:,1);
    
    WLF_all(ctr_s).rateA = tau_bA_grid(ctr_s);
    WLF_all(ctr_s).subsA = WlfA;
    WLF_all(ctr_s).subsB = WlfB;
    WLF_all(ctr_s).flag  = ind_neg_b;
    
end

fprintf(['     All    (1)Profit    Q_rdi    (2)RDi    Q_rde     (3)RDe   (4)Wdom   '...
        '(5)Wfrn   (6)Wint   (7)Trf   (1+2+3)  (1+2+3+5)  (1+2+3+4)  (1+2+3+6)\n\n']);

disp((WlfA./(Wlf_o*ones(1,size(WlfA,2)))).^(1/(1-psy)));
disp((WlfB./(Wlf_oB*ones(1,size(WlfB,2)))).^(1/(1-psy)));

CEQ = (WLFA./(Wlf_o*ones(1,size(WLFA,2)))).^(1/(1-psy));

[Val,Ind] = max(CEQ,[],2);
Optimal_rnd_over_time = [horizons' Ind tau_bA_grid(Ind)' (Val-1)*100];

% save Optimal_rnd_over_time_data CEQ WLF_all WLFA WLFB Wlf_o psy horizons tau_bA_grid MCUTM MCUTX

delete(gcp('nocreate'))
